Anomalous spatial diffusion and multifractality in optical lattices 
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Transport of cold atoms in sliallow optical lattices is characterized by slow, nonstationary momen- 
tum relaxation. We here develop a projector operator method able to derive in this case a generalized 
Smoluchowski equation for the position variable. We show that this explicitly non-Markovian equa- 
tion can be written as a systematic expansion involving higher-order derivatives. We use the latter 
to compute arbitrary moments of the spatial distribution and analyze their multifractal properties. 
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Slow relaxation is a hallmark of complexity. The 
absence of clear separations of time scales in complex 
systems often leads to nonexponential decay laws and 
anomalous diffusion properties Prominent systems 
exhibiting this kind of strange kinetics are cold atoms in 
dissipative optical lattices. Optical lattices are highly 
tunable periodic standing-wave potentials created by 
counter-propagating laser beams By decreasing the 
depth of the potential, the atomic dynamics can be mod- 
ified in a controlled manner from diffusive/ergodic to 
anomalous/nonergodic 0]. In the semiclassical regime, 
in which emission and absorption of single photons only 
slightly perturb the atomic motion, the latter can be de- 
scribed by a phase-space equation with nonlinear friction 
This Klein-Kramers equation has been successfully 
used in theoretical as well as in experimental inves- 

tigations 043- In particular, the predicted non-Gaussian 
momentum distribution Q, the divergence of its second 
moment in shallow potentials and the spatial su- 
perdiffusion [l^l have been experimentally observed. 
However, the general equation governing spatial diffu- 
sion, required to describe these experiments, is missing. 

In the theory of stochastic processes, the equation de- 
scribing the spatial motion is usually derived with the 
help of the standard procedure of adiabatic elimination of 
fast variables . This method exploits the rapid relax- 



ation of the momentum to a stationary state to extract an 
equation of the Fokker-Planck type for the slow position 
variable out of the full phase-space equation. The ob- 
tained equation is Markovian and does not depend on the 
initial condition. For shallow optical potentials, by con- 
trast, the momentum decay is algebraically slow [5| and 
the momentum distribution nonstationary [l^ — a well- 
defined relaxation time, hence, does not exist. As a re- 
sult, the standard elimination technique cannot be used. 
In the present paper, we develop a general projector op- 
erator method that is able to handle both slow momen- 
tum relaxation and nonstationarity. When applied to the 
optical lattice, we obtain a generalized non-Markovian 
Smoluchowski equation in the form of a systematic ex- 
pansion involving higher order derivatives. Similar ex- 
pansions, though Markovian and relying on stationarity. 



have been considered in the past |13l . Il4l | , in particular in 
the context of the Rayleigh model of Brownian motion 
with finite-time collisions 3, 1^ . They have been shown 
to provide more accurate descriptions than second-order 
equations [13,, jj| . Common to these studies is the pres- 
ence of a small parameter which allows the truncation of 
the expansion. Due to the inherent scale invariance of the 
atomic dynamics, such a small parameter does not exist 
in the case of the optical lattice. In the following, we use 
the generalized Smoluchowski equation to determine the 
generic moments of the spatial distribution, and reveal 
their multifractal structure. We further support our an- 
alytical predictions with detailed numerical simulations. 

Projector operator method. The Klein-Kramers equa- 
tion for the phase-space density W{x,p, t) is of the form 



dtW{x,p,t) = {L, + Lp)Wix,p,t), 



(1) 



where the linear operators and Lp are given by 

L,^-pd^ and Lp^dp{-F{p) + D{p)dp). (2) 

The semiclassical momentum-dependent drift and diffu- 
sion coefficients for atoms in the optical lattice are 0-01 , 



F{p) 



IP 



1 + [p/pcY 



D{p) = Do 



Di 



1 + {P/Pc 



(3) 



where the friction coefficient 7, the diffusion constants Dq 
and Di, and the capture momentum pc can be directly 
computed from the microscopic Hamiltonian of the prob- 
lem. Our goal is to derive the Smoluchowski equation for 
the position distribution, Wx{x,t) = J dpW{x,p,t). To 
this end, we introduce the projection operator P = 
such that Pf{p,x,t) = Wp{p,t*) J dp' f{x,p',t) for any 
integrable function f{x,p,t). Here Wp{p,t*) denotes 
the nonstationary solution of the momentum equation 
at time t* which we assume to be large, but fi- 

nite. In the quasistationary limit, r = < — i* <C i*, we 
have LpWp{p,t*) ~ 0. We check numerically the valid- 
ity of this approximation below. The function Wp{p,t*) 
is an even function of p for odd friction and even dif- 
fusion coefficients, like in Eq. Q, so that PL^P = 0. 
We define A{x,p,t) = PW{x,p,t) and B{x,p,t) = 



2 



(1 — P)W{x,p,t). The two functions A and B respec- 
tively belong to the subspace of relevant and irrelevant 
variables. The elimination method consists of two steps 
First, we apply the operators P and 1 — P to 
Eq. ([T|) to obtain two coupled equations for A and B. 
We then formally solve the equation for B and substi- 
tute the result in the equation for A, thus eliminating 
the irrelevant variables. This procedure is most conve- 
niently carried out in Laplace space. Using the notation 
f{x,p,s) = /g (iTe~'''^/(a;,p, r), we obtain the following 
equations for the Laplace transforms of A and B: 

sAis)-A{0) = PL,B{s), (4) 
sBis) - 5(0) = (L, + Lp)A{s) + [Lp + (1 - P)L,,]B{s). 

Solving the second equation for B{x,p, s), wc find 

sA{s)-A{0)^PL,H-^[L,Ais) + LpAis) + B{0)], (5) 

where we have defined H = s — Lp — [1 — P)Lx. Equation 
(0) is exact. For fast momentum decay, it is customary 
to approximate i/^^ ~ ~^p^ [Si- i^ST^e expand H^^ 
with respect to the coupHng term (1 — P)Lx |17| . 



H-' = j2[s-Lp]-(^^+'Hi-prL: 

and use the transformation, 

-1 n-l-l 



(6) 
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(7) 



The integral over r is finite, since stability requires Lp to 
have nonpositive eigenvalues. Further noting that (1 — 
P)" = (1 - P) (n > 0), we can write Eq. §1) as, 

Jo 

+ / dTT''eS^^-"'^'-{l- P)L^. (8) 

„=1 "'0 

For each value of n, we have to evaluate the three terms 
on the right-hand side of Eq. ([S]). The first term is 

PL^^^—^ J dTr"e(^^-"'^(l-P)L;'+iA(s) 



(9) 



where we have used the fact that e^''^p"^^Wp{p,t*) 
is the solution of the momentum equation with initial 
condition f{p) — p"+^W^p(p, i*). We have accordingly. 



= -d:+'W,ix,s)Wp{p,t*) J dTT^'e-'^ 
X / dp'p'''+^Wp{p'X)Ep{p',T), 



e'^-^p^+^WpipX) = Jdp'Wp{p,T\p',0)p'"+^Wp{p',t*) 
with the conditional probability density Wp{p,T\p' ,0). 
The second term vanishes for stationary processes. For 
nonstationary processes, it is of order r/t* , and can hence 
be neglected in the limit t <^t*. The third term can be 
set to zero by further assuming that the joint density fac- 
torizes at r = 0, W{x,p,Q) = Wp{p,t*)Wx{x,0). Since 
Ep{p,T) ~ J dp'p'Wp{p' ,t\p,0) is an odd function of p, 
Eq. ([HI is only non-zero for even n. Combining Eqs. 
^ and going back to time space, we eventually obtain 
the generalized, non-Markovian Smoluchowski equation, 

drW,Xx,T)^Yl / rfT'c„(r-T',r)9;w,(x,T'), (10) 

with the generalized correlation functions. 



c„(T,r) 



in -2) 



dpp''-^Ep{p,T)Wp{pX)- (11) 



For n = 2, Eq. (jlip gives the two-time momentum cor- 
relation function, C2(r, T) = {plr + t*)p(t*)). We stress 
that due to the nonstationarity of the momentum pro- 
cess, the correlation functions C„ depend explicitly on 
t* . Furthermore, Eq. (jTUl) reduces to the usual station- 
ary, Markovian, second-order Smoluchowski equation in 
the long-time limit, r ^ I/7, when the mean momen- 
tum relaxes exponentially fast, Ep{p,T) ~ e~^'^ . In this 
case, the spatial diffusion coefficient is = dr C2(r) 
ll|. For algebraically slow momentum relaxation with 
no characteristic decay time, the nonstationary, non- 
Markovian nature of Eq. ([TU)) will persist, even at long 
times, and the expansion cannot be truncated at the sec- 
ond order (see below). We can next derive an equation 
for the moments of the spatial distribution by multiply- 
ing Eq. (fTO|) by and then integrating over x. After 
integrating the nth term by parts n times, we find, 



9.(x™(t)) 

rn 

= E 



(m — n — 2)! Jq 



(12) 

dr'C„(r-T',r)(.T'"-"-2(r')) 



The mth moment is hence fully determined by the first 
(m — 1) moments and the correlation function (jlip . 

Application to optical lattices. Let us now apply the 
above formalism to the problem of transport of cold 
atoms in shallow dissipative optical lattices. For the drift 
and diffusion coefficients the asymptotic long-time 
behavior of the mean momentum Ep(p,t) is [l8| . 



Ep{p,-r) 



(4i:)oT)^"V"e 



2r(a + l) 



xM(|,a + l,^^ 



(13) 



where M{a,b,z) is the confluent hypergeometric func- 
tion, r(z) the Gamma function, and a = j/{2Dq) + 1/2. 
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For simplicity, we here set Di = and Pc = 1 (we 
will treat the general case elsewhere [l^). In addi- 
tion, the nonstationary, finite-time momentum distribu- 
tion Wp{p,t*) (the infinite covariant density) reads [l^ . 



Wpip,t*) 



(a>l) 

{a < 1) 
(14) 



with T{a, z) the lower incomplete Gamma fimction and 
Z = y^r{a - l)/T{a - 1/2) the partition function. We 
will use the above expressions to evaluate the asymptotic 
behavior of the spatial moments ([T^ to arbitrary order. 

The coefficients C„ (r, t* ) are directly related to the Tith 
moment of the position distribution via Eq. p2p . Evalu- 
ating the asymptotic behavior of Eq. ([TT|) with the help 
of Eqs. (fTSj) and we find the general expressions, 

(x"(T))-(x"(r))~C„,„ 

for a > n + 1 

3"+l_„ 



for 



1 < a < n + 1 



X < 



^*t+i-a^„ for 1< a < f + 1 



for a < 1 



Equation (|15p predicts normal behavior of the nth mo- 
ment for a > 71+1, anomalous, subballistic superdiffusion 
for n/2 + 1 < a < 71 + 1, and quasiballistic superdiffusion 
for a < n/2 + 1 and a < 1. The asymptotic formulas 
of the coefficients Cq,_„ can be computed analytically for 
all n using Eq. [Toj . We here only provide the two 
lower order terms. For n = 2, we obtain, 



.(4Do)--°r(. 2)r(3-o) ^4-. 2<«<3 

2Zr(a-i) r(5-a) 



(4^0 



Zr(a)(2-a) 

4Do(l -a)r r2 
while, for n = 4, we find 



for 1 < a< 2 (16) 
for a < 1 



{x\r)) - (xHn) 

2(4J>o)='-° 

zr(Q)(3-Q) 



^*3-a^4 



for 1 < a < 3 



(17) 



(4Do) (l-a)(2-a)r^r4 for a<l 



The second moment has been measured in the regime 
1 < a < 3 in Ref. and the transition from normal to 
superdiffusion has been observed. The case a < 1 has 
been examined recently, and ballistic behavior has been 
seen 0. In Figs, n and H we compare the theoretical 
results, Eqs. ([TBI) and (fTT)) . for the second and fourth mo- 
ment with Langcvin simulations of the full dynamics ([T]) 




FIG. 1. (color online) Second moment of the position distri- 
bution, Eq. H16p . for three different values of the parameter 
a (black solid). The colored dots are the results of Langevin 
simulations. Parameters are Do = 1 and t* = 5000. 
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FIG. 2. (color online) Fourth moment of the position distri- 
bution, Eq. (|17p . for three different values of the parameter 
a (black solid). The colored dots are the results of Langevin 
simulations. Parameters are Do = 1 and t* = 5000. 



(Euler method with time step 0.1 for 2.4 x 10^ trajecto- 
ries). The excellent agreement confirms the validity of 
the quasistationary aprroximation used in their deriva- 
tion. In both cases, the prefactor depends explicitly on 
the initial time t*, indicating the nonstationarity of the 
process. Let us emphasize two important points. First, 
the second moment is often calculated as a double time 
integral of the two-time momentum correlation function 
[lit . Trying to extend this approach to higher moments 
appears a formidable task since it involves the computa- 
tion of 7T,-time momentum correlation functions, contrary 
to the present method which only requires the general- 
ized two-time correlation functions C„. On the other 
hand, the nth moment of the equilibrium distribution, 
W^qip) = (l/2')(l-hp2)(i/2-"), diverges when n > 2a-2. 
Hence, had we used the latter in the definition of the pro- 
jector operator P, as is done in the standard elimination 
procedure [ll[, the corresponding correlation functions 
C„, Eq. pT|) . would not have been defined. 
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An essential property of the generalized Smoluchowski 
equation (jlOp is that at least n terms in the expansion 
are necessary to correctly compute the nth moment. For 
instance, terminating the expansion (|10p at the second 
order yields the fractional diffusion equation (2 < a < 3), 



(18) 



where we have introduced the Riemann-Liouvillc frac- 
tional derivative,^ ^r^/CO = ^/^i^) !ldt'{t-t'Y-^ f{t'), 
with < < 1 [20| • Equation (fT5|) describes a subballistic 
superdiffusive process. Its solution can be expresssed in 
terms of a Fox function and satisfies the scaling form (2lj 

W^[x,t) ^t^g{S,) with ^ = (19) 

The latter implies the scaling law for the moments: 

(|x|'?(t)) cx r^«. (20) 

Equation (j20p . however, only agrees with the general ex- 
pression (jl5p for the second moment. Higher-order mo- 
ments require higher-order terms to be taken into ac- 
count. These terms, in turn, will break the simple scaling 
of Eq. (fT5|) . leading to a complex multifractal structure. 

Multifractality. The concept of multifractality quan- 
tifies the variation of the scaling properties of a process 
[2^ . It often results from the overlapping of different 
scaling behaviors [2^ . Consider the gth moment , 



(|x|9(r)) cxT«(«). 



(21) 



A process is said to be monofractal, if the function ^(g) 
is linear, ^{q) = cq (normalization implies ^(0) = 0). For 
example, c = 1/2 for Gaussian diffusion, while c — 1 
for ballistic motion. For a multifractal process, ^(g) is a 
nonlinear function, indicating the presence of a contin- 
uous set of scaling exponents. Multifractality has been 
observed in many complex systems like, for example, tur- 
bulent and disordered mesoscopic systems (for a review 
see Ref. [2^). For integer moments, q = n, we can deter- 
mine the function £^{q) from Eq. pS]) . We find. 



for q < a — 1 



^ + l- Q;for a-l<q<2a-2 
q for g > 2q! — 2 or a < 1 . 



(22) 



Meanwhile, for nonintcgcr moments. Fig. [3] shows that 
^(g) is in general a nonlinear function of q. Hence, the 
exponent of the spatial gth moment exhibits a nontrivial 
dependence on q and a. For instance, for any given value 
of a, the moments (x^(r)) of order g < a — 1 behave 
as for normal diffusion. This is directly related to the 
finiteness of the equilibrium moment (p^'^). By contrast, 
for a— l<g<2a — 2, motion is subballistic superdiffu- 
sive, and even quasiballistic superdiffusive for q > 2a — 2. 
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FIG. 3. (color online) Exponent ^(g) of the moments 
(|a::|'(r)) oc r^^'^ as a function of the order q, for different val- 
ues of a. The solid line corresponds to normal Gaussian dif- 
fusion (5(g) = whereas the dash-dotted line represents 
quasiballistic expansion (£,{q) = q). The solid squares are 
obtained by fitting the long-time behavior of Langevin sim- 
ulations. The multifractality of the process is clearly visible, 
as 5(g) is not a linear function of g. It interpolates between 
normal diffusion for lower-order and quasiballistic motion for 
higher-order moments. Same parameters as in Fig. (1). 



The latter property illustrates the non-Gaussian nature 
of the underlying process: Even if the second moment 
increases linearly with time (a > 3), there always exist 
some higher-order moments which exhibit superdiffusive 
behavior. Thus, the second moment alone is not sufficient 
to fully capture the dynamics, as is typical for multifrac- 
tal processes. For a < 1, the process is monofractal. 

Conclusion. We have developed a general projector 
operator method to derive a spatial Smoluchowski equa- 
tion in situations where there is no time scale separation 
between position and momentum variables. Our formal- 
ism is able to treat both algebraic momentum relaxation 
and nonstationary momentum dynamics. We have shown 
that spatial diffusion is described by a nonstationary and 
non-Markovian generalization of the Smoluchowski equa- 
tion that can be written as a systematic expansion in- 
volving higher-order derivatives. We have applied this 
equation to the problem of transport of cold atoms in 
shallow dissipative optical lattices and obtained explicit 
expressions for spatial moments of arbitrary order. We 
have further established their generic multifractal behav- 
ior. Because of the latter, a truncation of the expansion 
is not possible. As a consequence, the lowest order term 
alone — given by a fractional diffusion equation — does not 
suffice to account for the complexity of the dynamics. 
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